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The  purpose  of  this  study,  the  results  of  which  are  reported  herein,  was 
to  formulate  a  mathematical  model  of  the  basic  VSTOL  jet,  and  to  validate 
its  appropriateness  by  performing  a  series  of  computational  experiments  on 
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finite  element  algorithm  for  solution  of  the  parabolized  form  of  the  time- 
averaged,  three-dimensional  Navier-Stokes  equations  (3DPNS).  The  essential 
ingredients  of  this  algorithmic  description  are  presented  and  discussed. 

A  production  finite  element  computer  code  (COMOC: 3DPNS)  was  utilized  to  conduct 
the  computational  experiments,  and  the  description  of  problem  definition  and 
data  deck  preparation  is  contained  at  the  end  of  the  report.  The  key  results 
of  the  computational  experiments  are  discussed,  regarding  the  basic  causal 
mechanisms  of  the  VSTOL  jet,  as  well  as  results  for  confirmation  test  cases 
documenting  viability  of  the  constructed  computational  model. 
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INTRODUCTION 


In  the  transition  from  hover  to  wing-borne  flight,  a  significant  portion  of 
VSTOL  aircraft  lift  is  furnished  by  engine  thrust.  The  injection  of  a  high- 
velocity  jet  of  exhaust  (or  air),  at  almost  right  angles  to  the  aerodynamic 
surface,  and  to  the  cross-flow  induced  by  the  aircraft  forward  flight,  pro¬ 
duces  an  extremely  complicated  three-dimensional  flowfield  which  can  signifi¬ 
cantly  impact  aerodynamic  performance.  Since  the  basic  phenomenon  of  trans¬ 
verse  jet  injection  is  so  fundamental  to  VSTOL  performance,  it  is  crucially 
important  that  a  firm  understanding  of  cause  and  effect  be  established.  The 
purpose  of  this  project  was  to  ascertain  dominant  fluid-dynamic  mechanisms 
characterizing  the  basic  VSTOL  jet,  using  a  computational  simulation  of  a  math¬ 
ematical  description  of  tiie  associated  three-dimensional  turbulent  flowfield. 

Since  the  transverse  jet  is  fundamental  to  many  physical  processes,  in 
addition  to  VSTOL  aircraft,  an  extensive  base  of  experimental  data  and  linear¬ 
ized  theoretical  analyses  is  in  existence.  The  first  consequential  theoretical 
study  was  reported  in  the  dissertation  of  Chang  [l] ,  who  employed  potential 
flow  theory  and  bound  vortex  filament  concepts  to  predict  the  shape  of  the 
separation  boundary  between  a  uniform  onset  flow  perpendicular  to  a  cylindrical 
jet.  A  series  solution  yielded  the  initial  roll-up  near  the  jet  exit  plane. 

A  (hand)  numerical  solution  produced  the  down-field  roll-up,  which  illustrated 
evolution  to  the  hallmark  horse-shoe  cross-sectional  shape.  The  basic  concept 
of  application  of  potential  flow  theory  has  been  extended  and  refined  to  a 
great  extent.  Viscous-corrected  potential  models  can  now  accurately  predict 
gross  far-field  flow  characteristics  such  as  jet  centerline  trajectory,  lateral 
spread,  cross-section  shape  and  mean  axial  velocity  at  each  cross-section, 
c.f.,[2,  3]. 

Jordinson  [4]conducted  pioneering  experiments  on  a  round  VSTOL  jet,  that 
confirmed  the  horse-shoe  cross-section  contours  predicted  by  the  potential  flow 
models.  However,  his  results  also  confirmed  that  the  VSTOL  jet  induced  the 
injection  plane  boundary  layer  flow  to  become  entrained  into  the  wake  region 
behind  the  jet.  This  action  was  not  predicted  by  the  elementary  theory,  which 
prompted  the  viscous-correction  procedures.  This  action  is  not  characteristic 
of  a  jet  issued  from  an  isolated  orifice  [5],  and  has  become  recognized  as  per¬ 
haps  the  consequential  contributor  to  alterations  of  VSTOL  aircraft  perfor¬ 
mance.  "Entrainment"  is  properly  interpreted  as  the  total  effect  of  vortex 
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roll-up  plus  turbulent  mixing,  while  "blockage"  describes  inviscid  flow  around 
an  equivalent  solid  body.  The  potential  flow  model  equivalent  is  therefore  a 
cylinder  with  suction  coupled  with  an  impirical  accounting  for  viscous  effects. 

These  observations  from  Jordinson's  data  prompted  conducting  of  a  wide 
range  of  experimental  tasks,  cf.,[6-12].  McMahon  and  Mosher  [8]  published 
photographs  of  oil  flow  streaklines  on  the  injection  plate,  see  Figure  1,  that 
provide  qualitative  evidence  of  the  difference  between  a  VSTOL  jet  and  the 
equivalent-diameter  solid  blockage.  For  the  latter,  Figure  la,  stagnation 
points  exist  both  upstream  and  downstream  on  the  center-line  and  are  connected 
by  an  apparent  inviscid  flow  streamline.  Outside  this  streamline,  the  exterior 
flow  simply  responds  to  the  blockage  induced  by  the  body.  The  flow  interior  to 
this  streamline  curves  inwards  to  the  centerline  and  appears  to  form  a  closed 
recirculation  zone.  Since  the  plate  is  impervious,  an  axial  flow  must  become 
induced,  directed  away  from  the  plate  and  parallel  to  the  solid  cylinder  for 
some  distance.  Thereafter,  the  shedding  of  the  Karman  vortex  street  probably 
results. 

The  streaklines  for  the  circular  VSTOL  jet,  Figure  lb,  are  consequentially 
different  (note  that  the  concentric  circle  is  part  of  the  experimental  appara¬ 
tus  and  not  a  streakline).  The  upstream  stagnation  region  does  not  appear  con¬ 
sequentially  different,  indicating  blockage-dominance.  However,  at  mid-jet, 
the  surface  flow  is  directed  almost  radially  inwards,  in  distinct  contrast  to 
solid  blockage  flow.  Two  streamlines  are  symmetrically  oriented  downstream, 
which  separate  the  complex  wake  interior  flow  from  the  deflected  free-stream. 
The  incoming  streamlines,  which  divide  the  flow  region  into  entrained  or  de¬ 
flected  segments,  intersect  the  downstream-dividing  streamlines  symmetrical ly, 
and  a  weak  stagnation  line  appears  to  connect  these  two  points  (essentially  on 
top  of  the  circle) . 

A  close-up  from  a  similar  test  [9]  is  shown  in  Figure  2,  which  clearly 
delineates  the  near  wake  streakline  patterns.  A  weak  stagnation  point  is 
indicated  about  one  jet  diameter  downstream  on  the  symmetry  plane,  with  a 
transverse  streamline  separating  wake  inflow  from  outflow.  The  dominant  ex¬ 
terior  flow  deflection  streamline  is  clearly  evident,  with  the  blockage  flow 
bifurcating  at  the  intersection  with  the  transverse  streamline. 
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Figure  1  Oil  Flow  Streaklines  On  The  Injection  Plate 
From  McMahon  And  Mosher  "8"  . 


Additional  insight  into  the  basic  transverse  jet  is  provided  by  experimen¬ 
tal  data  obtained  using  conventional  pressure  probes  in  the  farfield  [10-12]. 
Figure  3  illustrates  the  essential  VSTOL  geometry  and  coordinate  description, 
and  "far-field"  is  nominally  the  region  downstream  of  X^/D  *  5,  where  D  is  the 
jet  original  diameter,  flosher  [10]  concludes  that,  for  a  circular  jet  with 
velocity  ratio  4  £  A  £  10,  which  for  incompressible  flow  is  the  ratio  of  jet 

velocity  1).  to  the  onset  freestream  velocity  U  ,  entrainment  is  the  primary 
J  00 

bending  mechanism,  and  becomes  more  dominant  as  A  increases.  Kamotani  and 

Greber  [ll]  studied  the  A  =  8  circular  jet  and  pursued  considerable  data 
interpretation.  In  particular,  they  determined  the  locus  of  plane  orientations 
along  the  jet  path,  with  unit  normal  parallel  to  the  local  extremum  velocity, 
ie.  the  "jet  path."  In  Figure  3,  this  plane  is  spanned  by  the  X,,,  X^  coordi¬ 
nates,  and  the  normal  vector  is  parallel  to  X^.  Their  data,  plotted  in  the 
vertical  (YZ)  plane.  Figure  3,  produced  the  familiar  horse-shoe  profiles 
previously  determined  by  Jordinson. 

When  the  local  velocity  vector  is  resolved  into  scalar  components  in  the 
local  transverse  plane  coordinate  system,  the  jet  cross-section  shape,  and 
transverse  plane  counter-rotating  vortex  flow,  immediately  become  visible. 

For  example.  Figure  4  is  the  symmetric  half-plane  distribution  of  axial 
velocity  for  A  =  8,  at  stations  X^D  =  7  and  Xj/D  =  23.  The  characteristic 
"kidney"  shape  of  the  isovels  is  clearly  illustrated.  Since  the  jet  is 
essentially  axisymmetric  at  the  point  of  injection,  a  characteristic  action  of 
the  VSTOL  jet  appears  a  preferential  erosion  of  the  potential  core  in  the  near 
wake  region,  producing  a  double  maxima  off  the  symmetry  plane.  Figure  4c) 
shows  the  transverse  plane  velocity  vector  distribution  at  X^/D  =  23;  the 
center  of  the  vortex  is  coincident  with  the  extremum  axial  velocity,  Figure 
4b).  Transverse  velocity  data  for  the  A  =  8  circular  jet,  confirms  existence 
of  the  centered  double  vortex  structure  as  near  to  injection  as  X^/D  =  5.2, 

[12] .  These  data  have  been  examined  for  similitude;  Figure  5  summarizes  the 
minimal  data  spread  for  definition  of  jet  trajectory  over  a  range  of  velocity 
ratios. 

Limitations  with  conventional  probes  have  precluded  gathering  of  velocity 
data  closer  to  the  injection  plane  than  X^/D  =  5.  (The  emergence  of  the  three- 
color  laser  Doppler  velocimeter  will  perhaps  remove  this  obstacle.) 

Additional  characterization  of  the  basic  VSTOL  jet  has  thus  been  limited  to 
measurement  of  pressure  distribution  on  the  plate  forming  the  injection  plane 
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Figure  2  Oil  Flow  Streaklines  For  A  Circul 
From  Margason  And  Fearn  [9] . 


Figure  3  Coordinate  Description  For  VSTOL  Jet 


Experimental  Velocity  Distributions,  Circular 

Cross-Section  Jet,  X  =  8,  From  Kamotani  and  Greber  [l 


[10-14].  For  example.  Figure  6  compares  two  data  sets  for  the  injection  plate 
surface  pressure  coefficient  distribution,  for  the  circular  jet  at  X  =  8.  Only 
the  upstream  distribution  appears  analogous  to  potential  flow  about  a  solid 
cylinder.  The  extensive  low  pressure  region  downstream  of  mid-jet  is  a  dis¬ 
tinctive  characteristic  of  the  VSTOL  jet,  and  probably  a  (the)  significant  aero¬ 
dynamic  influence,  since  in  reality  the  injection  plane  is  an  aerodynamic  sur¬ 
face.  Alternative  jet  cross-sections  and  shapes  modify  the  details  of  the  asso¬ 
ciated  pressure  distribution,  see  Figure  7,  but  do  not  alter  the  basic  character. 

A  computational  simulation  of  the  causal  mechanism  of  the  extensive  low 
pressure  region  is  reported  [15].  For  this  analysis,  an  aerodynamic  contour  was 
assumed  to  exist  surrounding  the  jet,  and  separating  the  region  of  transverse 
flow  affected  primarily  by  blockage,  from  that  dominantly  influenced  by  entrain¬ 
ment.  This  contour  was  assumed  a  porous  Joukowski  airfoil,  with  the  surface 
onset  (suction)  velocity  distribution  an  input  parameter  of  the  simulation 
study.  A  series  of  computer  experiments  were  run  to  optimize  agreement  of  com¬ 
puted  pressure  coefficient  on  the  airfoil  surface,  with  experimental  data  on  the 
airfoil  image  on  the  injection  plate.  Figure  8  summarizes  the  results;  the  en¬ 
trainment  case.  Fig.  8b),  clearly  produces  a  large  aft  region  of  low  pressure, 
in  qualitative  agreement  with  experiment,  Fig.  6  and  7b).  In  comparison  to  the 
zero  entrainment  (solid  airfoil)  prediction.  Fig.  8a)  the  pressure  recovery  in 
the  trailing  edge  region  is  completely  absent. 

These  experimental  data,  and  the  computational  experiment,  appear  to  confirm 
that  the  large  low  pressure  region  is  the  direct  effect  of  the  entrainment 
action  of  the  basic  VSTOL  jet.  The  farfield  velocity  measurements  confirm  the 
preferential  wake  erosion,  of  the  initially  circular  jet  cross-section,  and  the 
evolution  of  a  transverse  vortex  pair.  The  axial  velocity  contours  flatten 
broadside  to  the  onset  flow  to  produce  a  bluff  rather  than  aerodynamic  cross- 
section.  While  these  data  provide  a  rather  comprehensive  characterization,  the 
dominant  causal  influence  remains  no  more  well  defined  than  being  an  interaction 
between  the  jet  and  the  cross-flow. 

The  purpose  of  this  study,  the  results  of  which  are  reported  herein,  was  to 
formulate  a  mathematical  model  of  the  basic  VSTOL  jet,  and  to  validate  its 
appropriateness  by  performing  a  series  of  computational  experiments  on  the 
discrete  analog  (numerical)  approximation  to  the  mathematical  description.  Since 
the  ideal  VSTOL  jet  problem  is  essentially  incompressible,  steady,  turbulent  and 
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Figure  6  Experimentally  Measured  Surface  Pressure  Coefficient 
Distributions  On  Injection  Plane,  Circular  Jet,  X 


fully  three-dimensional,  the  mathematical  description  must  be  quite  comprehen¬ 
sive.  The  approach  selected  was  to  utilize  a  continuity-constrained,  finite 
element  algorithm  for  the  parabolized  form  of  the  time-averaged,  three-dimen¬ 
sional  Navier-Stokes  equations  (3DPNS).  The  essential  ingredients  of  this 
algorithmic  description  are  discussed  in  the  section  on  Problem  Definition.  A 
production  finite  element  computer  code  (COMOC: 3DPNS)  was  utilized  to  conduct 
the  computational  experiments,  and  the  description  of  problem  definition  and 
data  deck  preparation  is  contained  at  the  end  of  this  report.  The  key  results 
of  the  computational  experiments  are  discussed,  regarding  the  basic  causal 
mechanisms  of  the  VSTOL  jet,  as  well  as  confirmation  test  case  results  docu¬ 
menting  viability  of  the  constructed  computational  model. 


a)  Bluff  Rectangular 


c)  Streamline  Rectangular 

Figure  7  Experimentally  Measured  Surface  Pressure  Coefficient 

Distributions  On  Injection  Plane,  Various  Jets,  \  =  4,  From  ; 14] 


/ 


/ 

a)  Solid  Blockage  With  Zero  Entrainment 


b)  Circular  Jet  With  Entrainment,  A  =  8 

Figure  8  Computed  Pressure  Coefficient  Distribution  At  Injection  Plane, 
From  Baker,  Manhardt  And  Yen  [15  1 
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PROBLEM  DEFINITION 


Overview 

The  requirement  is  to  establish  a  partial  differential  equation  system  pro¬ 
viding  an  adequate  description  of  the  basic  VSTOL  jet  configuration.  As  men¬ 
tioned  in  the  Introduction,  the  essential  features  are  three-dimensional  and  tur¬ 
bulent,  with  the  flow  basically  steady  and  of  constant  density.  The  complete 
three-dimensional,  time-averaged  Navier-Stokes  equations  (3DNS)  provide  the 
accurate  description;  however,  their  numerical  solution  is  essentially  intrac- 
tible  with  present  computer  hardware.  The  three-dimensional  parabolic  Navier- 
Stokes  equations  (3DPNS)  is  a  simplification  of  the  3DNS  system,  that  is  appro¬ 
priate  for  steady  flow  prediction  in  aerodynamic  configurations  wherein  the  pre¬ 
dominant  velocity  component  does  not  suffer  reversal,  and  where  certain  other 
requirements  are  met.  The  next  level  of  simplication  to  3DNS  reduces  the 
problem  description  to  two-dimensional,  which  is  unacceptable.  The  3DPNS 
system  is  numerically  tractible,  and  therefore  has  been  employed  for  the  com¬ 
putational  analysis  of  the  VSTOL  jet  problem. 

This  Section  presents  the  construction  of  the  3DPNS  equation  system  for 
the  VSTOL  jet  problem,  including  a  turbulence  closure  model  for  the  Reynolds 
stress  tensor.  The  finite  element  numerical  solution  algorithm  for  3DPNS  is 
described,  including  the  differential  constraint  procedure  employed  to  enforce 
the  (non-parabolic)  continuity  equation.  Documentary  numerical  results  are 
briefly  outlined  that  validate  the  theoretical,  accuracy  and  convergence  aspects 
of  the  resultant  computational  simulation  capability. 

Three-Dimensional  Parabolic  Navier-Stokes  Equations 

The  3DPNS  equation  set  is  a  simplification  of  the  steady,  three-dimensional 
time-averaged  Navier-Stokes  equations.  In  Cartesian  tensor  notation,  and  em¬ 
ploying  superscript  tilde  and  bar  to  denote  mass-weighted  and  conventional  time¬ 
averaging,  respectively  [16] ,  the  conservative  equation  form  for  an  isoenergetic 
fluid  is 


1(3 }  ■tar'jj =  0 


(i) 


L ( ou . )  3  4 — 1 3u •  u  ■  +  5i  ■  .  +  ru-'u'  -  r  •  ■ 

1  3xi L  1  j  i  j  ij_ 

J 


=  o 


(2) 
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L(3k)  « 


*  iUiUj  ixj 


T  P£  =  0 


L(St)  *  kh‘  *  ! +  ci  w  C 


♦c?r"0  <«: 

In  equations  1-2,  p  is  (constant)  density,  u.  is  the  mean  velocity  vector,  p  i; 

J 

pressure,  5..  is  the  Kronecker  delta.  The  Stokes  stress  tensor  c.  ■  is 

J  ’  J 

defined  in  terms  of  the  Reynolds  number  Re  as 

°ij  =  0v  ^ij  ’  I  eij*W/Re  (5; 

and  puTuT  is  the  Reynolds  stress  tensor.  The  fluid  kinematic  viscosity  is 

^  v 

v,  and  E--  is  the  mean  flow  strain  rate  tensor 

•  J 

l  =  2Mi  +  iii 

ij  _  3Xj  3x.j  (6 

Equations  3-4  are  the  transport  equations  for  turbulent  kinetic  energy 
and  isotropic  dissipation  function,  as  obtained  using  the  closure  model  of 
Launder,  Reece  and  Rodi  [17]  for  the  pressure-strain  and  triple  correlations, 
and 

k =  itipir 


-  |  3Uj  3U-i  . 

■ 5  3  b  ;xkj:jk 


The  various  coefficients  are  model  constants,  Hanjalic  and  Launder  [18] 

P 

The  parabolic  Navier-Stokes  equation  set  is  derived  from  equations 
1-4  assuming  the  ratio  of  extremum  transverse  mean  velocity  component  to 
downstream  (axial)  component  is  less  than  unity,  and  that: 
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1.  the  downstream  velocity  component  suffers  no  reversal, 

2.  diffusive  transport  processes  in  the  downstream  direction  are 
higher-order,  hence  negligible,  and 

3.  the  overall  elliptic  character  of  the  parent  three-dimensional 
Navier-Stokes  equation  can  become  enforced  through  construction 
of  a  suitable  pressure  field. 


Viewing  Figure  3,  for  the  VSTOL  jet  problem,  the  x.\  (curvilinear)  coordinate 
defines  the  predominant  mean  flow  direction,  with  scalar  velocity  component  uj 
of  order  unity,  i.e.,  0(1).  Hence,  assume  0(u2)  a  0(5)  «-  0(u3)  and  0(5)  <  0(1). 
Then,  the  continuity  equation  1  asserts  that  the  downstream  variation  in  u3 
must  be  of  order  equal  to  appropriate  transverse  plane  variations  of  u2  and  u3; 


hence,  for 


3x  i 


*  h, 


Determination  of  the  relative  order  of  terms  in  the  momentum  equations  2 


is  straightforward.  For  the  Ui  equation,  since  O(pu^u^)  must  be  0(6),  the 

3  _ _  1  J 


The  assumption  that  the  xj 


term  (piTfuT)  is  higher  order  and  discarded 

diffusion  is  negligible  permits  setting  t-  (En)  =  0;  further,  -r- 

f3u2'1  3X1  9X2 


3__ 
3x  3 


3x 


iJ 


3u2 

1>Xi 


-  0(6) 


Deletion  of  these  terms  is  the  fundamental  step  to  the  parabolic 


approximation,  as  it  removes  the  elliptic  character  in  the  downstream  direction. 
The  convection  term  (puiUi)  instills  the  initial  value  form  that  permits 

d  X  i 


marching  the  solution  for  u3  in  the  downstream  direction, 
form,  denoted  Lp(.),  for  the  Ui  momentum  equation  is 


The  final  3DPNS 


lp(^.)  *jz,Gwi)  +-ff1+  T5T2  Luru£ 


u 3 !  =  0 

J  (9) 

As  a  final  note,  should  x.  correspond  to  a  curvilinear  coordinate  description, 

J 

the  derivatives  expressed  in  equation  9  become  the  covariant  derivative. 

In  agreement  with  the  parabolic  concept,  the  order  of  pressure  variation 
in  the  transverse  plane  must  be  determined  by  the  lowest  order  terms  appearing 
in  equation  2  written  on  u2  and  u3.  Each  transverse  derivative  of  pu2u^  and 
pu3u;j  is  0(1),  while  all  other  terms  are  0(6)  and  higher. 

V 
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The  resultant  three-dimensional  parabolic  approximation  is  made  numerically 
tractible  by  taking  the  divergence  of  the  approximation  to  equation  2.  Retain¬ 
ing  the  higher-order  convection  and  laminar  diffusion  terms  for  generality,  the 
3DPNS  equivalent  of  both  transverse  momentum  equations  is  [19], 


L("p)  =  iS  - 


y  iL 


:  U  .  U  •  +  o  U'  u  ■ 

L  j  c  J 


(10) 


which  defines  the  3DPNS  limited  index  summation  convention,  1  £  (i,j)  £  3, 

2  <  l  <_  3. 

Equation  10  defines  an  elliptic  boundary  value  problem  for  determination 
of  pressure  distributions  in  the  transverse  plane.  The  pressure  field  that 
satisfies  this  Poisson  equation  consists  of  complementary  and  particular 
sol utions . 


P(Xj)  *  pc(*()  ♦  fpUj)  (u) 

The  comp! ementary  solution  is  assumed  to  satisfy  the  convection  pnenomena. 


,2n  U  ) 

1  n  1  =  +  - i-J- 

3X2  3X.3X 

K  J  ~ 


=  0 


(12) 


subject  to  the  Dirichlet  boundary  conditions  known  for  p(xi,x£).  Elsewhere, 
an  appropriate  boundary  condition  for  pc  is  homogeneous  Neumann.  The  particular 
pressure  pp  satisfies  equation  10,  less  the  convection  term,  subject  to  homo¬ 
geneous  Dirichlet  boundary  conditions  on  boundary  segments  where  pc  is  known. 

The  critical  aspect  affecting  application  of  the  3DPNS  equation  set  to 
analysis  of  the  VSTOL  jet,  is  knowledge  of  the  boundary  values  for  p  (xi,  x^), 
as  required  via  assumption  3  of  the  3DPNS  argument.  Viewing  Figure  3,  the 
VSTOL  jet  is  a  fully  three-dimensional  problem  with  elliptic  coupling  through¬ 
out  in  its  entirety.  However,  see  Figure  5,  the  trajectory  of  the  jet  for  the 
first  few  jet  diameters  away  from  the  injection  plane  is  nominally  vertical. 
Hence,  if  the  ‘computational  box"  surrounding  the  jet  and  the  near  field  flow 
is  sufficiently  large  in  lateral  (X,Y)  extent,  it  is  fair  to  assume  that  the 
potential  flow  pressure  exists  on  the  box  boundaries,  and  for  some  vertical 
distance  (parallel  to  the  Z  axis).  Therefore,  the  validity  of  the  3DPNS 
analysis  is  expected  to  be  limited  to  a  region  close  to  the  injection  plate 
surface,  eg.,  0  <  Xi/D  <  4,  provided  the  lateral  boundaries  are  sufficiently 
remote.  For  the  analyses  reported  herein,  these  boundaries  of  the  3DPNS 
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simulation  were  placed  at  (x  |/D  -  3.  The  jet  path  was  assumed  straight  and 
perpendicular  to  the  plane.  The  complementary  pressure  solution  boundary 
condition  p,  on  the  upstream  and  lateral  side  boundaries,  were  set  by  the 
farfield  potential  solution  for  flow  about  a  cylinder.  The  complementary 
pressure  pc  at  the  outflow  boundary  was  determined  using  Bernoulli's  equation, 
which  admitted  reversal  of  the  transverse  flow  and  was  compatible  with  the 
lateral  farfield  potential  solution. 


Reynolds  Stress  Closure  Model  for  3DPNS 

A  closure  model  for  the  kinematic  Reynolds  stress  -uTuT,  appearing  in 
equations  3-4  is  required  to  complete  the  3DPNS  order  of  magnitude  analysis. 
A  stress-strain  rate  constitutive  equation.  Baker,  et.  al.  [20]  establishes 
the  lead  terms  of  the  kinematic  form  as 


-Vj  =  -^ij 


.is § 


j 


C;cXe  E  ; 

=.  I  K  Kj 


(13) 


E.  ■  remains  the  symmetric  mean  flow  strain-rate  tensor  given  in  equation  6. 

^  3 

This  expansion  results  from  re-expression  of  triple  correlations  within  the 
Reynolds  stress  transport  equation  using  the  model  of  Launder,  Reece  and  Rodi 
[I7j,  and  is  a  general ization  of  the  original  analysis  by  Gessner  and  Emery 
[21].  In  equation  13,  a-jj  is  a  diagonal  tensor  in  principal  coordinates. 


(uk“k>V;u  <14) 

The  a.j  are  coefficients  admitting  anisotropy,  where  a3  =  Ci,  and  a2  =  C3 
=  a3.  The  Ca  are  defined  as,  see  Launder  et.  al.  [17]. 

2 2(Coi  -  1)  -  6(4Cq;  -  5) 

Cl  3  33 (Co ~  -  2C0 z ) 

4(3Ca2  -  1) 

°2  H  11 (Co  1  -  2Cc : ) 

22 ( C o  i  -  1)  -  1 2 ( 3C o 2  -  1) 

Cj  5  33 (Co i  -  2C o ; 7 

44C,:  *  22CoiCo;  -  128C0;  -  36Co:  +  10 
C-  S  1 6 5 ( Co  l  -  2Co2 ): 
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where  C0i  and  Ca2  are  "universal"  empirical  constants.  Suggested  values  are 
Co i  3  2.8  and  Z02  -  0.45,  Hanjalic  and  Launder  [18] .  In  rectangular  Cartesian 
coordinates,  and  retaining  terms  of  the  first  two  orders  of  significance,  the 
kinematic  Reynolds  stress  scalar  components  are 

0(6)  0(6:) 


—y—?  ~  .  r  r  k'  ]  ■  SUl  “  .  I  jUI  •  -  ; 

-ufuf  *  Cik  -  C:C4—  .-r~\  +  •-rr-\ 

£  J  X  2  J  ljX3j  J 


-UjU2*  =  C3k  -  C;Cu— |j 
-W-  =  C3k  -  C2C4|r;| 

— r  k:;r3Ul 


r  ,  2 


-U1U3  = 


k3  f3Ux] 
e‘  L3X2J 

2  V  2  f'  i'i  1 

l'  :  0  li  2  i 

-  2C4— 

-  L3X2- 

k  3  foil 

12  V  2  1*,'  .1 

„  or  -  1 ~u  -  . 

£-  i_3 X 3 J 

t  ;_dX3^ 

!~3u  i~ 

+ 

3  U  ;  1 

L3X2. 

-  4£“ !  3X5  ;.-'X3 

~-X2j 

+  -fo1  r 3Ui  +  3u2) 

3X2 (5Xi  3X2J 

ilk! 

-  c.c44^^  + 

’bx3j 

‘  3X21.5X3 

3X:J 

-u£uj  = 


rr  k3 ! 5Ui  ;U 1 
*  ~  |_3 X 2  5  X  3_i 


+  o3U: '  5Ui  + 

5X 3  V3X j.  aX3j 

r  k2  1  3U;  ,  5U3  1 

L“  £  '  SX3  3X2 : 


With  equation  16,  the  ordering  of  terms  in  equations  3-4  can  be  completed 
to  establish  the  appropriate  30PNS  approximation  as,  1  <_  i  <^3,  2  <_  l  <  3, 


l/OO  -^<50,10  ♦  ^ 


+  c  + 


+  =  0 


l^u)  *trcsfufu: 


k  — r— -  3e 


+  cl  pUTU;  +  c2  3  f  =0 
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Finite  Element  Solution  Algorithm 

The  consistently  ordered  3DPNS  equation  system  has  been  identified.  For 
the  dependent  variable  set  q j ( )  e  (q }  =  (us,  u2,  u3,  p,k,e}T,  the  governing 
system  includes  equations  1,  9,  10,  13,  17  and  18.  Equations  9,  10,  17  and  18 
exhibit  the  initial  value  term  that  permits  the  space-marching  of  3DPNS. 
Equation  13  exhibits  elliptic  boundary  value  character  with  parametric  initial- 
value  dependence.  The  continuity  equation  1  solution  will  become  recast  and 
utilized  as  a  differential  constraint. 

The  generalized  form  of  the  30PNS  description  is 

*  k(=°‘V  *  kkqj  *  fuj  *°  us 

where  f . .  and  s.  are  specified  non-linear  functions  of  their  arguments,  as 

X.  J  J 

determined  by  the  index  j.  The  solution  domain  Q  is  defined  as  the  product  of 
R2  and  Xi,  for  all  elements  of  xi  belonging  to  the  open  interval  measured  from 
xi(0),  i.e., 

n  =  R2  x  Xi  =  {(x.,x:):  x^eR2  and  Xie[xj (0) ,x:) 

The  boundary  9S1  is  the  product  of  the  boundary  9R  of  R2  and  Xi,  9Q  =  9R  x  Xj. 
Thereupon,  the  generalized  differential  boundary  constraint  is 


*ai  =  0 


where  the  a^  are  specified  coefficients  and  fL  is  the  outwards  pointing  unit 
normal  vector.  Finally,  an  initial  distribution  for  qj  on  f20  =  R2  x  xi(0) 
is  required. 


VVXl)  i  q:(xj 


The  finite  element  numerical  solution  algorithm  for  equations  19-22 
defines  the  approximation  qj(x^,xi),  to  the  (unknown)  exact  solution  qj(x^,xi), 
which  is  constructed  from  members  of  a  finite  dimensional  subspace  of  Ho(fO> 
the  Hilbert  space  of  all  functions  possessing  square  integrable  first  deriva¬ 
tives  and  satisfying  equation  21.  Hence, 


9j(x^,Xi)  ~  Qj(x^,Xi)  =  s_-j  qjlx^xO 
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and  the  elemental  approximation  is  defined  as 


dj (x.  jX].)  =  {Nk(x,)}T{QJ(x1)}e  ^24 

L. 

In  equations  23-24,  j(J)  is  a  free  index  denotir^  members  of  { qn } ,  and  sub- 
or  super-script  e  denotes  pertaining  to  the  e^-  finite  element,  Si  =  x  X[. 

y  6  G 

The  elements  of  the  row  matrix  ^(x^)}  are  (linear,  k  =  1)  polynomials  on  x^, 
2  <  £  <  3,  constructed  to  form  a  cardinal  basis  [22]. 

L. 

The  functional  algorithm  requirement  is  to  render  the  error  in  minimum 
is  some  norm.  This  is  accomplished  using  the  finite  element  algorithm,  by  re- 

L. 

quiring  the  generated  errors  Lp(q^)  and  £(q.)  to  be  orthogonal  to  the  function 

h  J  J  D  -h 

space  employed  to  define  q.  ,  and  that  the  discrete  approximation  L/(p  )  to 

J 

the  continuity  equation  1  be  enforced  as  a  differential  constraint.  Identify¬ 
ing  the  multiplier  set  8..,  these  constraints  are  combined  [22]  to  form  the 
theoretical  statement  of  the  finite  element  solution  algorithm  as 


R2{Nk}LP(qJ)dx+  8:/3R{Nk}£(qjh)dx  + l-/R27{Nk}LP(3h)dx  =  {0} 


(25) 


Equation  25  defines  a  system  of  ordinary  differential  equations  written 
on  the  jet  direction  x2 . 


[C]{QJK  +  [U]{QJ>  +  C  FLO  ]  (QL }  +  (SJ'=  {0}  (26) 

Using  the  trapezoidal  integration  rule,  and  substituting  equation  25  yields. 


iFJ}  s  -  {QJ}j  -  xEQ0>j+1  +  (QJ)j]  =  -0) 


(27) 


which  defines  a  system  of  non-linear  algebraic  equations  for  determination  of 
the  elements  of  { Q J ( x i ) } .  The  Newton  iteration  algorithm  for  equation  27  is 


[J(FJ)]P+1(cQJ}P:;  -  *(FJ}P+1 
written  on  the  iteration  vector  {<5QI},  where 


(28) 


w;-;:’,  -=  *  «j># 


(29) 


The  remaining  issue  is  the  non-parabolic  continuity  equation  1,  which 
governs  first-order  effects  on  the  mean  velocity  field  Q. .  Since  the  formal 
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3DPNS  momentum  equation  9  is  written  on  Ui  only,  =  { u2 » u 3  >  is  required  to  be 
determined  from  the  solution  of  equation  1.  This  is  accomplished  in  the  finite 
element  algorithm  formulation  by  enforcing  a  measure  of  the  solution  of  equation 
1  as  a  differential  constraint  on  the  solution  of  the  (higher-order)  3DPNS  approx¬ 
imation  to  equation  2  written  on  u^.  Retaining  the  highest  two  orders  of  terms, 
the  3DPNS  order  of  magnitude  analysis  for  transverse  momentum  equations  yields 


+ 


ix,L-u?.uk 


(30) 


Note  that  equation  30  is  of  the  form  of  equation  19,  and  employs  the  3DPNS 
limited  index  2  <_  k  j<  3.  The  middle  two  terms  in  the  second  bracket  are  0(1), 
while  the  remaining  terms  are  all  0(6).  The  boundary  condition  statement  21 
is  appropriate,  upon  the  retention  of  terms  of  0(62)  in  the  Reynolds  stress 
equation.  The  initial  condition  statement  for  equations  30  is  expressed  in 
equation  22. 

The  finite  element  algorithm  statement  for  equation  30  is  given  by  equation 
25,  upon  specification  of  the  form  of  the  term  modified  by  $2.  The  theoreti¬ 
cal  concept,  borrowed  from  the  variational  calculus  [22],  is  to  enforce  a 
measure  of  the  continuity  equation  (solution)  as  a  differential  constraint  on 
solution  of  the  transverse  momentum  equations.  This  solution  measure  must  span 
R2,  and  must  vanish  as  continuity  becomes  satisfied;  the  appropriate  measure  is 
the  harmonic  function  4>(x^),  the  solution  to  the  Poisson  equation 


'  2. 


ip(3)  =  -  frr(3Qi )  E  0 


(31) 


The  boundary  conditions  for  <j>  are  homogeneous  Dirichlet  everywhere  at  farfield. 
Equation  31  becomes  homogeneous  as  the  continuity  equation  1  becomes  satisfied, 
and  <j>  becomes  null  as  a  consequence  of  the  boundary  condition  specifications. 

A  detailed  discussion  of  the  algorithmic  embodiment  of  the  differential  con¬ 
straint  concept  is  given  in  [22], 
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Documentary  Results 

The  finite  element,  differential-constraint  numerical  solution  algorithm 
for  the  3DPNS  equation  system  1,  9,  10,  13,  17,  18,  30  and  31,  for  a  turbulent, 
subsonic  three-dimensional  analysis  in  bounded,  semi -bounded  and. fully  unbounded 
solution  domains  ft  ,  is  operational  in  the  C0M0C:3DPNS  computer  program.  A 
latter  report  Section  presents  the  detailed  data  deck  set-up  procedure  for  the 
VSTOL  jet  simulation  study.  Additional  documentation  on  data  preparation  for 
C0M0C.-3DPNS  will  become  available  [23]. 

Documentation  of  theoretical  concepts,  and  numerical  tests  to  evaluate 
accuracy  and  convergence  are  reported  elsewhere,  cf.  [19],  [22].  As  becomes 
apparent  in  the  results  discussion,  the  essential  requirements  for  the  3DPNS 
VSTOL  simulation  are  robust  enforcement  of  (the  non-parabolic)  continuity 
equation,  and  accounting  of  the  turbulence  interaction  between  the  jet  and  the 
cross-flow. 

Test  case  results  attest  to  excellent  performance  of  the  algorithm  in  C0M0C: 
3DPNS  in  each  area.  The  two-dimensional  form  of  the  3DPNS  algorithm,  ie. ,  2DPNS, 
must  accurately  predict  ducted  and  boundary  layer  flows,  as  well  as  free  shear 
layers.  Table  1,  from  [19],  summarizes  the  comparison  between  the  direct 
boundary  layer  solution,  and  the  continuity-constraint  2DPNS  algorithm  results, 
for  a  laminar,  incompressible  zero-pressure  gradient  boundary  layer  flow.  The 
digit  of  significance  ranges  over  six  digits  (106)  and  the  agreement  is  excel¬ 
lent.  The  2DPNS  solution  is  actually  more  accurate  near  the  wall  (x2/S  = 

0.0009),  where  a  better  approximation  to  vanishing  x2-  derivative  is  achieved. 

The  intrinsic  global  measure  of  error  in  exact  conservation  of  mass  is  the 
energy  norm, 

E(*V)  =  |^-  |~  d?  (32) 

Rn 

where  4>h  is  the  approximate  solution  to  equation  31.  The  evaluation  of  E ( - ,  • ) 
decreases  monotonical ly  during  iteration  at  j+1,  and  for  this  laminar  flow 

“8  “5 

check  case,  E(*,*)  £  10  for  algorithm  convergence  at  e  =  10  .  The  corres¬ 

ponding  levels  for  a  turbulent  boundary  layer  or  ducted  flow  are  E ( * , * )  <  10 
for  e  =  10"“  [19]. 
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The  critical  evaluation  of  the  Reynolds  stress  closure  system,  equations 
13-16,  is  also  reported  in  [19].  For  three-dimensional,  incompressible  flow  in 
a  straight  channel  of  uniform  rectangular  cross-section,  an  axial  vortex  pair 
is  induced  in  each  corner  for  a  turbulent  flow,  while  none  is  generated  if  the 
flow  is  laminar.  The  results  of  the  3DPNS  simulation  confirm  these  assertions 
on  both  counts.  In  particular,  generation  of  the  counter-rotating  vortex  pair 
for  turbulent  flow  was  computationally  traced  to  non-isotropy  of  the  transverse 
plane  normal  stresses  -u'u'  .  Viewing  equation  16,  this  influence  is  produced 
by  a  term  of  0(62),  which  confirms  the  importance  of  a  consistent  ordering  pro¬ 
cedure.  Figure  9  shows  the  excellent  qualitative  comparison  of  the  transverse 
plane  velocity  vector  distribution  as  produced  by  C0M0C:3DPNS  [19],  and  the 
experimental  data  of  Melling  and  Whitelaw  [24]. 

Table  1 

Transverse  Velocity  Distributions,  u2(x2)xl03 
Laminar  Incompressible  Boundary  Layer 


Coordinate 

_(x2/<5) 

Boundary 

Layer  Solution 

Continuity 
Constraint  Solution 

0.0 

.0 

.0 

0.0009 

.0000011 

.0000001 

0.0021 

.0000054 

.0000030 

0.0035 

.0000149 

.0000114 

0.0095 

.00011 

.00011 

0.031 

.00118 

.00118 

0.10 

.0128 

.0127 

0.67 

.139 

.138 

1.0 

.218 

.218 
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RESULTS  AND  DISCUSSION 


Overview 

The  requirement  of  the  3DPNS  computational  simulation  of  the  VSTOL  jet  is 
to  assess  fluid  dynamic  factors  dominating  formation  of  the  counter-rotating 
axial  vortex  apir,  see  Figure  6,  and  entrainment  of  the  cross-flow  into  the 
jet  wake  region,  cf.  Figures  2-3.  The  3DPNS  solution  domain  was  defined  as  the 
symmetric  half-plane,  with  the  circular  jet  located  mid-domain.  Figure  10a. 

For  reference.  Figure  10b)  is  the  transverse  plane  potential  velocity  vector 
distribution  used  to  initialize  u^(xi(o),  x^)  at  the  nodes  of  ft(o).  The  tail 
of  each  vector  is  at  a  node,  the  non-uniform  discretization  is  constituted  of 
M  =  576  triangular  finite  element  domains  R2.  Using  this  discretization,  C0M0C: 
3DPNS  required  approximately  200,000  words  of  central  memory  to  solve  for  eight 
degrees  of  freedom/node,  using  the  described  implicit  algorithm,  equations  26- 
29.  No  exterior  memory  was  utilized. 

Boundary  and  Initial  Conditions 

The  computational  simulation  requires  solution  of  a  non-linear  system  of 
ordinary  differential  equations.  Hence,  the  members  Ui,  u^,  k  and  c,  must  be 
defined  initially  at  the  nodes  of  $1(0).  Secondly,  since  the  finite  element 
algorithm  statement  has  transformed  the  elliptic  boundary  value  character  of 
the  3DPNS  equations,  boundary  conditions  on  all  members  of  (QI(xi)}  are  re¬ 
quired  specified  everywhere  on  3ft  =  3R2  x  Xi.  Table  2  summarizes  the  boundary 
condition  statement  for  each  variable,  on  segments  A-D  of  the  boundary  3R  , 
see  Figure  10a).  Basically,  line  AB  is  a  symmetry  plane,  CDA  is  the  farfield 
potential  boundary,  and  BC  is  inflow/outflow  with  vanishing  normal  derivatives. 

The  specification  of  suitable  initial-conditions  {QI ( 0 ) >  is  a  perplexing 
problem,  since  no  experimental  data  are  available  for  guidance.  For  the  simu¬ 
lations  reported,  the  initial  jet  was  assumed  to  be  either  of  circular  or 
square  cross-section,  with  a  constant  axial  velocity  (distribution)  ui(x?)  =  U^. 
The  initial  ut  distribution  outside  the  initial  jet  cross-section  was  assumed 
constant  at  a  level  smaller  than  the  initial  jet  velocity  U..  This  constant 

J 

was  selected,  dependent  upon  the  particular  test  and  the  associated  numerical 
stability.  The  constant  is  required  non-zero,  to  prevent  the  3DPNS  equation 
set  from  becoming  singular,  dependent  upon  the  imposed  cross-flow  velocity 
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Table  2 


Boundary  Condition  Specifications  for 
3DPNS  Simulation  of  VSTOL  Jet  in  Cross-Flow 


Variable  Equation  Boundary  Segment  Boundary  Condition 


ABCDA 


Zero  normal  derivative 


BC 

CDA 

Qa  30  ABC 

CDA 

k  17  ABC 

CDA 

e  18  ABC 

CDA 

Pc  12  ABC 

CDA 

Pp  10  ABC 

CDA 

<P  31  AB 

BCDA 


Vanishing  normal  derivative 
Specified  cross-flow 
Vanishing  normal  derivative 
Specified  cross-flow 
Vanishing  normal  derivative 
Fixed  at  freestream  level 
Vanishing  normal  derivative 
Fixed  at  freestream  level 
Vanishing  derivative 
Fixed  at  potential  cross-flow 
Vanishing  derivative 
Zero 

Vanishing  derivative 
Zero 
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level  (cf.,  the  order  of  magnitude  analysis).  The  VSTOL  jet  velocity  ratio 

is  thereby  X-  U./U  . 

j  m 

For  tests  conducted  with  zero  cross-flow,  the  initial  distribution  for  the 
transverse  plane  velocity  field  (x i »x^ )  was  zero.  With  >  0,  and  following 

considerable  numerical  experimentation,  the  potential  cylinder  flow  field  shown 
in  Figure  10b)  was  used  to  define  0^(0)  outside  the  initial  jet.  Inside  the 
initial  jet,  u^(0)  5  0.  Thus,  Bernoulli's  equation  yields  the  corresponding 
initial  distribution  of  pc  on  £2(0). 

Establishing  suitable  initial  distributions  for  turbulent  kinetic  energy 
and  dissipation  function  is  more  perplexing  since  the  interaction  on  the  up¬ 
stream  face  of  the  jet  is  extremely  complicated.  Since  the  project  scope  is 
to  assess  dominant  factors,  and  in  the  total  absence  of  data,  a  step  distribu¬ 
tion  in  turbulence  levels  was  assumed  appropriate.  The  initial  level  of  tur¬ 
bulent  kinetic  energy  k  was  assumed  a  distinct  constant  inside  and  outside  the 
jet  over  £2(0).  The  initial  levels  of  dissipation  were  specified  unique  con¬ 
stants,  as  well,  to  produce  distinct  levels  of  turbulent  kinematic  viscosity 
=  Ci*k2/e  inside  and  outside  the  jet,  eg.,  f-  . 

A  series  of  computational  experiments  were  required  to  ascertain  numerically 
acceptable  levels  of  k(0)  and  e(0),  within  the  constraint  1  <  v^/v^  <  102.  The 
k  and  e  equations  17-18  are  directly  coupled  in  the  non-homogeneous  annihilation 
terms  pc  and  Pe2/k,  and  the  farfield  solutions  depend  strictly  upon  this  balance 
(since  all  spatial  gradients  essentially  vanish).  Boundary  levels  for  k  and  c 
were  fixed  on  segments  CDA,  Figure  10a),  and  the  3DPNS  algorithm  executed  to 
determine  solutions  of  equations  17-18  in  the  farfield.  The  new  levels  were 
then  set  on  CDA,  and  the  tests  repeated  until  the  interior  solution  and  boundary 
data  agreed.  For  =  10,  the  suitable  levels  were  determined  as  k(0)  =  0.0005 
and  e(0)  =  0.0012.  For  example,  then  for  *  102,  k(0)  =  0.005  and  e(0)  = 

0.012  are  acceptable  initialization  levels  inside  the  jet. 

Validation  Numerical  Results 

As  discussed  in  the  following  sub-section,  the  3DPNS  simulation  of  the 
VSTOL  jet  produces  a  transverse  plane  flowfield  in  substantial  qualitative 
agreement  with  anticipated  results  and  extensions  of  the  modest  data  base. 

The  requirement  is  to  ascertain  that  the  simulation  is  free  from  numerical  dis-  j 

turbances  that  would  by  themselves  produce  such  results.  For  this  purpose,  j 

i 

I 

I 

i 
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several  test  cases  were  executed  in  the  absence  of  an  imposed  crossflow,  ie. 

=  0,  to  quantify  algorithm  performance  with  regards  to  symmetries  and  robust 
conservation  of  mass.  The  latter  is  particularly  critical  since  the  sole  mech¬ 
anism  for  generation  of  u^  f  0  is  through  continuity  which  is  employed  as  the 
differential  constraint  on  eouation  30  that  otherwise  exhibits  a  null  solution. 

Examples  were  conducted  for  initially  circular  and/or  square  cross-section 
jets.  Figure  11  shows  the  3DPNS  computed  transverse  velocity  distributions  on 
0.5  £  Xj/D  £  2.0,  for  a  jet  of  initially  square  cross-section,  for  =  0  and 
L).:U.  =  1.0:0. 2,  where  U.  is  the  Ui(0)  inside  the  jet  and  LIT  is  Ui(0)  outside 

J  J  J  J 

the  jet.  The  flow  Reynolds  number,  based  upon  jet  velocity  (37  m/s)  and  jet 
diameter  (width)  is  Re  =  .6  X  10",  and  the  turbulent  kinematic  viscosity  v*"  » 
102  was  held  a  constant.  The  computed  velocity  field  exhibits  essentially 
exact  symmetry.  Strong  flow  opposition  is  computed  on  the  boundary  of  the  jet 
at  xi/D  =  0.5,  and  thereafter  becomes  more  distributed  throughout  the  initial 
"potential  core."  Note  also  the  lateral  spreading  of  the  jet  boundary,  as 
defined  by  the  opposing  velocity  vector  locus.  All  velocity  vector  plots  are 
scaled  on  the  maximum  scalar  component  of  u^(xi,  x£).  The  caption  u™  under 
each  figure  legend  is  the  ratio  of  this  maximum  to  U..  Hence,  note  the  ex- 

_  m  J 

tremum  u^  was  reached  at  x^D  =  1.0,  u^  =  0.026,  and  thereafter  decays  by  30% 
in  reaching  xx/D  =  2.0. 

Figure  12  shows  the  3DPNS  computed  u^  distributions  on  0.5  £  x^D  £  1.0,  for 
the  case  in  Figure  11  but  with  U  - ; U  T  =  1.0:0.02,  ie.,  the  farfield  background 

J  J 

for  ui  outside  the  jet  is  1/50  of  Uj .  Comparing  to  Figure  11,  a  significantly 
stronger  entrainment  field  is  induced  at  Xi/D  =  0.5  (u™  =  0.086),  but  there¬ 
after  there  is  no  consequential  difference.  This  test  addresses  the  issue  of 
the  3DPNS  order  of  magnitude  analysis,  wherein  0(6)  <  0(1)  by  some  unspecified 
amount.  For  the  test  of  Figure  11,  0(6)  =  0(Uj/u™)” 1  K  0( 10" 1 ) ,  while  for 
Figure  12,  0(6)  =  0(0.02/0.05)" 1  =»  0  (1).  For  the  cases  with  non-zero  cross- 

flow,  where  U.:U  =  1.0:0. 1,  these  results  indicate  that  for  the  farfield  level 
j  °° 

of  U-  *  U^,  the  3DPNS  algorithm  should  function,  although  in  this  region  0(5)  * 
0(1)  is  not  robust  adherence  to  the  ordering  analysis.  Some  detailed  solution 

differences  would  occur  for  Uj  not  a  uniform  constant;  the  evolution  of 
( x x )  would  remain  essentially  unchanged,  however,  based  upon  these  results. 

The  final  validation  is  quantization  of  the  algorithm  capability  to  predict 
turbulent  shear  layer  interaction  on  a  grid  discretization  as  coarse  as  shown 
in  Figure  10a).  For  this  purpose,  the  two-dimensional  algorithm  (2DPNS)  was 
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c)  Xj/D  =  1.5 

uj  =  0.021 
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d)  xj/D  =  2.0 
uj  =  0.017 


Figure  11.  3DPNS  Computed  Transverse  Plane  Velocity  Vector  Distributions, 
Square  Jet,  U./UT/U^  =  1.0/0. 2/0.0,  /  =  102 . 
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b)  Xi/D  =  1.0 
u™  =  0.052 


Figure  12.  3DPNS  Computed  Transverse  Plane  Velocity  Vector  Distributions 
Square  Jet,  U./lT/U^  =  1.0/0.02/0.0,  vl  =  102 . 
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executed  to  simulate  a  turbulent  slot  jet  flow  at  U.:if  =  1.0:0.02,  using  a 

J  J 

uniform  discretization  of  R1  of  mesh  measure  equal  to  that  for  the  smallest 
element  domains  shown  in  Figure  10a),  as  located  on  the  boundary  of  the  jet. 
Figure  13  summarizes  the  2DPNS  solution  on  0  <  xj/D  £  1.0,  for  a  unit  Reynolds 
number/Length  of  Re  =  0.2  X  107/m.  The  action  of  turbulent  viscosity  smoothly 
erodes  the  slot  jet  boundary.  Entrainment  from  the  farfield,  u2  <  0,  and 
definition  of  the  jet  boundary  by  the  sign  change  in  u2,  is  clearly  illustrated 
in  Figure  13b).  The  initial  distributions  k(0)  and  e(0)  were  determined  from 
elementary  mixing  length  theory  [16] .  The  solutions  for  k  and  e.  Figures  13c)- 
d)  show  the  rapid  growth  in  both  variables  at  the  edge  of  the  slot  jet.  The 
resulting  decay  of  these  peak  levels,  as  the  strain  rate  su!/3x2  diminishes 
(hence  also  the  non-homogeneous  source  term  in  the  k  and  e  equations  17-18),  is 
also  evident. 


Figure  14  shows  the  corresponding  computed  distributions  of  the  Reynolds 
stress  tensor,  equations  13-16.  The  axial  normal  stress  ufuf  exhibits  larger 
gradients  and  peak  levels,  in  comparison  to  the  transverse  normal  stress  u2u2 , 
due  both  to  the  difference  between  Cx  and  C3,  and  the  action  of  the  0(52)  term. 
This  latter  influence  is  documented  [20]  of  dominant  importance  in  promoting 
prediction  agreement  with  detailed  experimental  data  for  an  airfoil  trailing 
edge  wake  flow.  The  Reynolds  shear  stress  and  turbulent  kinematic  viscosity 
distributions  exhibit  similar  trends. 


Circular  VST0L  Jet  Simulation  Results 


The  principal  subject  is  analysis  of  a  circular  cross-section,  subsonic 

VST0L  jet,  issued  perpendicular  to  a  flat  plate  into  a  subsonic  crossflow  at 

velocity  ratio  A  =  U./U  =  10.  The  results  of  the  validation  test  cases,  and 

J  00 

additional  numerical  experiments,  provided  the  initialization  procedure.  The 

initial  plateau  distributions  of  axial  velocity  u^O)  and  the  reference  cross- 

flow,  were  set  at  U.:U]:U  =  1.0:0. 2:0.1,  with  U-  =  36.6  m/s  (120f/s).  In  the 
J  J  00  J  _ 

absence  of  better  guidance,  the  initial  velocity  U.  exterior  to  the  jet  was 

J 

assumed  a  uniform  constant,  to  preclude  occurrence  of  fictitious  farfield 
gradients.  The  initial  distribution  for  transverse  velocity  1^(0)  is  potential 
flow  about  a  cylinder.  Figure  10b).  The  initial  levels  for  turbulence  kinetic 
energy  k(0)  ^  0.0005  and  dissipation  e(0)  <_  0.0012  were  set  at  distinct  constant 
levels  inside/outside  the  initial  jet  to  produce  the  desired  initial  levels  of 
1  <  vt  <  100.  The  reference  Reynolds  number  based  on  jet  orifice  diameter  is 
Re  =  0.6  X  104. 


30 


a)  Axial  Velocity  Ui 


b 


c 


d 


e 


b)  Transverse  Velocity  u2 


c)  Turbulent  Kinetic  Energy  k 


d)  Dissipation  Function  e 


Figure  13.  2DPNS  Computed  Solution  Field,  Turbulent  Rectangular 

Slot  Jet  Flow,  Uj/Uj  =  1.0/0.02,  xi/D  =  (0,  .25,  .5,  .75,  1.0}  = 

{a,  b,  c,  d,  eh 
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A 


a)  Normal  Stress  uTuT 


b)  Normal  Stress  ulul 


c)  Shear  Stress  uTuT 


d) 


Turbulent  Viscosity  v 


t 


Figure  14.  20PNS  Computed  Reynolds  Stress  Distributions, 

Turbulent  Rectangular  Slot  Jet  Flow,  Xi/D  =  {0,  .25,  .5,  .75,  1.0}  = 


{a,  b,  c,  d,  e}. 
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The  necessary  comparison  results  were  for  this  specification  in  the  absence 
of  cross-flow,  ie.,  Uro  =  0.  Figure  15  summarizes  the  3DPNS  solution  for  trans¬ 
verse  velocity  distributions  on  0.25  £  Xi/D  _<  1.25  for  laminar  flow.  Figure  16 
shows  the  comparison  turbulent  flow  results  for  vt(0)  *  102  inside  the  jet  and 
v^O)  *  10  everywhere  exterior  to  the  jet.  Both  figures  illustrate  the  contin¬ 
uity-induced  generation  and  entrainment  of  a  farfield  transverse  velocity  dis¬ 
tribution.  The  extremum  velocity  magnitude  occurs  on  the  jet  boundary  for  both 
cases;  however,  the  laminar  farfield  magnitude  is  considerably  smaller,  in  com¬ 
parison  to  its  extremum,  than  for  the  turbulent  case,  see  Figures  15c)  and  16c). 
Of  greater  fluid  dynamic  significance  regarding  entrainment,  the  extremum  tur¬ 
bulent  transverse  velocities  are  a  factor  of  20  or  more  larger  than  the  lami¬ 
nar  flow  extrema.  This  ratio  is  a  quantitative  indicator  of  the  enhanced  en¬ 
trainment  efficiency  of  the  turbulent  jet.  Figure  17  shows  the  3DPNS  solution 
distributions  for  =  {ui,  k,  e}  at  xi/D  =  1.25,  for  the  turbulent  circular 
jet  in  zero  cross-flow.  The  contour  levels  are  labeled  and  shown  in  the  appro¬ 
priate  legend.  The  distributions  for  Ui  and  k  are  close  approximations  to 
axisymmetric,  with  a  slight  flattening  indicated  along  the  45°  and  135°  rays 
from  the  center.  This  is  a  grid  induced  effect,  primarily  the  consequence  of 
the  initial  condition  approximation  on  the  union  of  rectangles,  Figure  10a). 

This  can  be  observed  in  the  smallest  contour  level  for  e  =  0.020,  Figure  17c), 
which  lies  just  above  the  background  initialization  level  of  e(0)  =  0.012. 

Since  the  application  of  +  0  will  destroy  any  point  symmetry  anyway,  these 
results  are  taken  as  confirmation  that  the  selected  combination  of  discretiza¬ 
tion,  boundary  constraints  and  initial  conditions  for  the  3DPNS  VST0L  jet  simu¬ 
lation  will  not  induce  a  spurious  flowfield  skewing  of  measurable  consequence. 

Figure  18  shows  the  3DPNS  algorithm  computed  evolution  of  the  transverse 

velocity  field  on  0  <  xi/D  <  1.25,  for  the  circular,  turbulent  VST0L  jet  in  a 

cross-flow,  U.:UT:U  =  1.0:0. 2:0.1.  The  extremum  transverse  scalar  component, 

J  J  00 

upon  which  each  figure  is  scaled  is  noted  in  each  legend.  The  blockage  effect 
of  the  jet  is  immediately  evident.  Figure  18b),  as  is  the  beginning  of  entrain¬ 
ment  in  the  lateral  farfield.  The  transverse  velocity  field  begins  to  penetrate 
the  downstream  interface  of  the  jet  by  xi/D  =  0.5,  Figure  18c),  and  a  decrease 
in  the  outflow  in  the  wake  region  is  evident.  This  becomes  considerably  pro¬ 
nounced  by  xi/D  =  0.'75,  Figure  18d),  where  wake  flow  reversal  has  begun. 
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Figure  15.  3DPNS  Computed  Transverse  Velocity  Distributions, 
Laminar  Circular  Jet,  Zero  Cross-Flow. 
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Figure  16.  3DPNS  Computed  Transverse  Velocity  Distributions, 
Turbulent  Circular  Jet,  Zero  Cross-Flow. 
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a)  Axial  Velocity  ut 


{a,  b,  d,  d,  e}  = 

{.3,  .5,  .7,  .9,  1.0} 


e  d  c  b  a 


b)  Turbulent  Energy  k 

{a,  b,  c,  d}  = 

{.013,  .025,  .040,  .044} 


i 

abcdcba 


c)  Dissipation  Function  <c 

{a,  b,  c,  d}  ~~ 

{.02,  .2,  .5,  .65} 


abcdc  ba 

Figure  17.  30PNS  Computed  Solution  {Q I > , 

Turbulent  Circular  Jet,  Zero  Cross-Flow,  xt/D  =  1.25 
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The  entrainment  from  the  lateral  farfield  has  also  become  stronger  with  pene¬ 
tration  on  the  lee  side  of  the  jet.  The  vortex  structure  interior  to  the  in¬ 
itial  jet  becomes  fully  developed  by  xi/D  =  1.0,  Figure  18e).  The  combined 
actions  of  blockage,  entrainment,  wake  flow  reversal  and  axial  vortex  are  each 
matured  and  clearly  evident  in  the  solution  at  Xi/D  =  1.25,  Figure  18f). 

The  results  of  Figure  18  can  be  interpreted  as  a  consequential  indication 
that  the  3DPNS  simulation  procedure  exhibits  potential  for  prediction  of  the 
near-field  essential  characteristic  action  of  the  elementary  VSTOL  jet  in 
cross-flow.  In  particular,  the  region  of  reversed  flow  and  angle  of  the  wake 
streamline.  Figure  18f),  are  in  qualitative  agreement  with  the  oilflow  streak¬ 
line  experiment  results.  Figures  1-2.  The  generation  of  the  axial  vortex 
(pair)  is  in  qualitative  agreement  with  farfield  data  [l0-12]see  Figure  4.  The 
jet  boundary  is  not  impervious  to  the  cross-flow  and  the  intrusion  of  entrain¬ 
ment  is  predicted  on  the  downstream  face  of  the  jet.  However,  as  clearly  em- 
phasi2ed  in  the  preceeding  discussions,  these  results  are  consequentially 
influenced  by  the  many  decisions  and  compromises  required  to  complete  the  math¬ 
ematical  specification,  in  particular  the  initial  conditions  and  size  and  re¬ 
finement  of  the  computational  solution  domain.  These  are  detailed  aspects  that 
require  individual  attention.  The  impetus  to  attack  these  subjects  is  hope¬ 
fully  enhanced  by  the  encourageing  results  of  this  prediction. 

Additional  computational  tests  were  conducted  to  quantize  the  influence  of 
gross  turbulence  factor  modifications.  Figure  19  summarizes  the  results  for 
the  VSTOL  jet  specification  of  Figure  18.  Figures  19a ) -b )  show  the  transverse 
velocity  field  u^  on  0.5  <_  Xi/D  £  1.0  computed  holding  v*"  s  10,  a  constant 
throughout  the  entire  solution  domain.  Some  farfield  entrainment  action  occurs, 
but  the  nearfield  crossflow  appears  almost  negligibly  deflected.  Figure  19c) 
shows  the  comparison  solution  at  x^D  =  0.5,  for  constant  v*  but  with  h  100. 
The  jet  is  considerably  more  impervious  to  the  crossflow.  In  comparison  to 
Figure  18c),  a  somewhat  greater  deflection  of  the  upstream  farfield  has 
occurred,  and  downstream  penetration  is  essentially  absent.  For  the  variable 
turbulence  simulation,  Figure  18,  an  extremum  vt  ~  160  was  computed  on  the  jet 
boundary  at  Xi/D  =  0.5.  Clearly,  the  turbulence  phenomena  exerts  a  predomi¬ 
nant  influence  on  predicted  results  which  adds  further  emphasis  to  the  need  to 
obtain  appropriate  high  quality  experimental  data. 
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Figure  19  3DPNS  Assessment  Of  Gross  Turbulence  Modifications, 

Circular  VST0L  Jet,  U.UT/U  =  1.0/0. 2/0.1 
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Figure  20  shows  the  3DPNS  solution  distributions  for  q?  =  {uj,k,  e}  at 
x,/D  =  1.25,  for  the  turbulent  circular  jet  with  cross-flow,  (J.-:UT:U  =  1.0: 
0.2:0. 1.  (The  transverse  velocity  solutions  are  shown  in  Figure  18.)  The  con 
tour  levels  are  labeled  and  numbered  in  the  appropriate  legend.  The  exact  com 
parison  solution  for  zero  cross-flow  is  shown  in  Figure  17,  and  distinct  dif¬ 
ferences  are  noted.  In  general,  the  essential  axisymmetry  for  zero  cross-flow 
has  been  altered,  with  cross-flow  applied,  and  for  xi/D  =  1.25,  the  profile 
transverse  spans  are  quite  flattened  with  the  isoclines  somewhat  distended  up¬ 
stream  along  a  nominal  135°  ray  from  the  circle  center.  The  extremum  level 
contours  for  k  and  e  both  occupy  a  larger  region  of  R2,  generally  within  the 
interval  135°  <  a  <  180°.  These  results  appear  the  combined  influence  of  the 
lateral  constraint  influence  of  cross-flow,  counter-balanced  with  the  effects 
of  entrainment  and  core  vortex  convection  to  distend  the  flow  profiles  prefer¬ 
entially  in  the  upstream  quadrant. 


a)  Axial  Velocity  uj 
{a,  b,  c,  d,  e)  = 
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Turbulent  Energy  k 
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c)  Dissipation  Function  e 
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abCdCba 


Figure  20.  3CPNS  Computed  Solution  { QI } 

Turbulent  Circular  Jet  With  Cross  Flow,  X  =  10, 
Xj/D  =  1.25. 
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SUMMARY  AND  CONCLUSIONS 


1 


A  mathematical  model  has  been  constructed  for  analysis  of  the  near-field 
evolution  of  the  basic  VSTOL  jet  in  an  applied  subsonic  crossflow.  The  resul¬ 
tant  numerical  approximation  has  utilized  a  continuity-constrained,  finite 
element  algorithm  for  solution  of  the  parabolized  form  of  the  time-averaged, 
steady  three-dimensional  Navier-Stokes  equations  for  turbulent  flow.  A  Reynolds 
stress  constitutive  equation  was  employed  for  turbulence  closure  which  required 
numerical  solution  of  the  parabolized  three-dimensional  form  of  the  turbulent 
kinetic  energy  and  dissipation  function  transport  equations.  Various  documen¬ 
tor  tests  were  executed  to  validate  correct  operation  of  the  C0M0C:3DPNS  com¬ 
puter  program  embodiment  of  the  theory. 

The  results  of  the  computational  simulations  of  the  VSTOL  jet,  with  and 
without  application  of  a  subsonic  crossflow,  have  yielded  results  of  technical 
interest  and  anticipated  merit.  Of  greatest  significance,  the  constructed 
3DPNS  model,  when  operated  without  fictitious  constraints  on  the  turbulence 
field  (evolution),  produced  a  flowfield  solution  for  the  circular  jet  in  sub¬ 
stantial  qualitative  agreement  with  the  available  sparse  data.  Without  initial- 
condition  generated  bias,  and  using  the  most  elementary  solution  starting  con¬ 
ditions,  the  3DPNS  model  predicted  the  essential  evolution  character  in  sub¬ 
stantial  completeness.  In  particular,  the  simulation  predicted  lateral  en¬ 
trainment,  axial  vortex-pair  initiation  and  inducement  of  the  wake  flow  into 
the  jet  region,  antiparallel  to  the  initialized  potential  crossflow  direction. 
These  documented  features  of  the  VSTOL  jet  were  generally  lost  when  the  tur¬ 
bulence  field  was  artificially  constrained,  indicating  that  the  characteristic 
action  is  a  turbulence-dominated  effect. 

It  should  be  emphasized  as  well  that  none  of  the  flow  characteristics  were 
consequentially  captured  without  robust  enforcement  of  the  (non-parabol ic)  con¬ 
tinuity  equation.  The  developed  constraint  algorithm  met  the  detailed  mathe¬ 
matical  requirements  and  accurately  enforced  the  detailed  solution  aspects 
related  to  the  3DPNS  ordering  analysis.  In  particular,  the  finite  element 
based  algorithm  maintained  the  (energy  norm)  error  in  exact  satisfaction  of 
continuity  at  E(.,.)  <_6.  X  10"5  for  algorithm  convergence  set  at  e  =  10  . 

The  3DPNS  algorithm  averaged  four  iterations  per  step  for  convergence  following 
a  few  extra  iterations  to  homogenize  initial  condition  error.  The  solution  on 
0  <_  Xj/D  <  1.25  required  approximately  60  CPU  minutes  to  execute  on  a  CYBER  175 
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computer,  using  200K  single  precision  words  of  central  memory. 

These  results  are  interpreted  as  highly  encouraging  with  respect  to  consti¬ 
tuting  a  basic  proof  of  concept.  Clearly,  the  utilized  discretization  is  too 
coarse  by  an  order  of  magnitude,  and  the  transverse  domain  span  too  small  by  a 
factor  of  two  or  more.  The  total  absence  of  any  data  to  help  generate  initial 
condition  distributions  is  certainly  a  liability.  Nevertheless,  the  numerical 
simulation  capability  is  three-dimensional  and  turbulent,  and  has  provided 
some  valuable  and  interesting  insight  to  quantify  dominant  aspects.  The  dis¬ 
cretization  and  domain  span  constraints  are  strictly  computer  details  (that 
can  be  solved  by  utilizing  the  CYBER  20X,  for  example).  The  3DPNS  algorithm 
can  be  directly  extended  to  compressible  and  non-isoenergetic  flows,  hence  be 
rendered  applicable  to  the  more  practical  problem  involving  heated  jets  (ex¬ 
haust).  The  Reynolds  stress  closure  model  could  be  replaced  by  solution  of 
the  parabolized  Reynolds  stress  transport  equations,  or  even  more  elegant 
procedures,  if  and  when  it  is  deemed  necessary  and/or  appropriate.  Each  of 
these  issues  constitutes  a  particular  detail,  building  upon  the  basic  structure 
of  the  computational  model  and  computer  program.  The  3DPNS  model  could  be  of 
considerable  usefulness,  for  engineering  analysis,  pending  emergence  of  the 
complete  3DNS  numerical  solution  capability  (and  the  CYBER  20X2  or  NASF). 
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APPENDIX 
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DATA  DECK  SPECIFICATIONS 


Input  facilities  for  the  COMOC  IV  Parabolic  Navier  Stokes  computer  program 
are  sophisticated  and  greatly  simplify  data  deck  preparation  and  modification. 
The  program  sequentially  scans  the  data  deck  and  operates  on  command  data  cards 
as  they  are  encountered.  Numerical  data  required  for  each  command  operation  are 
input  in  free  format  on  cards  directly  following  the  command  card.  Command  oper¬ 
ations  can  cause  arrays  to  be  filled,  initiate  a  series  of  solution  operations 
or  specify  output  formats  and  titles.  Command  card  sequence  is  quite  flexible 
and  care  has  been  taken  to  ensure  that  most  operations  which  must  be  performed 
sequentially  are  specifiable  under  one  command  name. 


Most  numerical  data  may  be  input  in  free  format.  Data  delimiters  may  be 
blanks  or  commas,  thus  allowing  for  esthetic  and  meaningful  arrangement  of  data 
and  simple  addition  or  deletion  of  interspersed  numbers.  Several  features  of 
free  format  input  which  greatly  simplify  repetitive  and  sequential  data  specifi¬ 


cations  are: 


Repetitive  Numbers: 
Fills  Array 

Repetitive  Sequence 
"(One  per  card  only) 

Fills  Array 

Skip  P  locations 
Fills  Array 

Increment  by  a  constant 
Fills  Array 

Exponential  Notation 
Fills  Array 


12.5*7 

12.  7.  7.  7.  7.  7. 
2(5.  2.  4. 


5.  2.  4.  5.  2.  4. 

10.  12.  3*P  22.  T 
10.  12.  V  V  V  22. 

5*50  10  T 
10  60  110  160  210 

6.  10.0E-2  14.0E-4 

6.  0.1.0014 


The  data  deck  for  a  VST0L  jet  in  a  crossflow  solution  is  segmented  into 
several  sections  for  description,  but  appears  in  sequential  order  as  indicated 
by  line  numbers.  The  data  has  been  lumped  into  meaningful  categories,  each  of 
which  can  be  related  to  specific  differential  equation  solution  components, 
(e.g.,  initial  condtions,  boundary  conditions,  etc.).  The  Fortran  MAIN  program 
contains  global  dimensioning  and  initialization  of  the  primary  differential 
equation  coefficient  arrays  and  parameter  lists.  A  listing  of  MAIN  together 
with  subroutines  GETTPR  and  N0DPPR,  which  provide  flow  initial  conditions 
specific  to  the  jet  crossflow  problem,  follow  the  data  deck  description. 


I.  NAMELIST  (Integer  Scalar)  Input 


230PNS 
3FKNAHE 
4  SNAME01 


5 

NODE  =»  330. 

6 

NM  -  3* 

7 

NEC)  =  9. 

8 

NEQAV2  =  If 

9 

NPVSX  =  30f 

10 

IBLAS  -  2. 

11 

NIMPLT  =  lr 

12 

KNTPAS  =  4  f 

13 

SEND 

Line 

Conmand 

2 

3DPNS 

3 

FEBANE 

4 

&NAME01 

5 

NODE 

LCOL 

KROW 

6 

NDP 

7 

NEQ 

NEQKNN 

NEQADD 

8 

NEQAV2,3 

NCNADD 

9 

NPVSX 

ITKE 

10 

IBLAS 

IU0NLY 

11 

NIMPLT 

NBAND 

NCNTIT 

12 

KNTPASS 

NPRNT 

NTPRNT 

NC 

13 

&END 

T  READ  NAMELIST  PARAMETER  DATA  -  INTEGER .REAL 


LCOL  =  30 f 

NDP  =»  10. 

KROW  = 

30. 

NEOKNN  =  5. 

NEQADD 

=  -4. 

NEQAV3  =  1, 

NCNADD 

=  5. 

NTABpr  =  3. 
IUONLY  =  3. 

ITKE  = 

1. 

NBAND  =  29. 

NCNTIT 

=  0. 

NPRNT  =  132. 

NTPRNT 

=  99999 

NC  =  8  > 


Solution  procedure  (3D  Parabolic  Navier-Stokes) 

Command  to  initiate  Namelist  read 

Integar  Namelist  Data 

Slightly  larger  than  the  number  of  nodes 
in  the  solution 

Larger  than  the  number  of  nodes  along  abcissa 

Larger  than  the  number  of  nodes  along  the  ordinate 

No.  of  Differential  Equations  plus  No.  of 
parameter  arrays  required 

No.  of  Dependent  Variables 

No.  of  Differential  Equations  solved 

No.  of  Differential  Equations  initially  not  solved 

Set  to  1  to  integrate  variables  U2  and  U3 

Begin  U3  integration  after  NCNADD  steps 

Number  of  points  in  the  input  Pc  Table 

Flag  indicating  solution  of  the  turbulent 
shear  stresses 

U2»  U3  convergence  multipliers  in  EPS 

Number  of  steps  until  convergence  on  UI  only 

Use  implicit  Newton  Iteration  procedure 

Maximum  bandwidth  of  Jacobian  Matrix 

Pass  counter  to  start  perturbation  pressure  solution 

Maximum  number  of  integration  steps  between  prints 

Number  of  columns  on  a  line  of  output 

Exclude  integral  parameter  print 

Output  format  number  field  width 

End  of  Integer  Namelist  Data 
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NAMELIST  (Real-Scalar)  Input 


14 

&NAME02 

15 

UINF  3  120.. 

TOFINF  3  533.. 

REFL  3  .0328084. 

16 

GUMULT  =1., 

17 

OSG  3  1 . . 

PPFACT  3  1., 

PCFACT  3  1., 

18 

VLDMLT  =  1 . . 

VCMULT  3  1.. 

RARK'AY  <  385  >  3  1.1 

IV 

RHSCAL  =  0. . 

U2STRS  3  1.. 

RHQIM  3  1.0. 

20 

C4EDSW  3  .03. 

OSUSQ  3  l.E-4. 

21 

SIMPLT  =l.E-5. 

TMULT  3  1.21. 

EFMULT  3  .01. 

22 

TO  =  1.0E-5. 

TD  3  10.0. 

DELP  3  20.. 

23 

CHIEPS  =1 ,E-4 . 

HMAX  3  4.0. 

24 

SEND 

1 .E-10, 


ALC  =  1.0. 
HSINIT  =  l.E-4, 


Line 

Command 

14 

&NAME02 

Real  Namelist  Data 

15 

UINF 

Velocity  non-dimensionalizing  factor 

TOFINF 

Temperature  non-dimensionalizing  factor 

REFL 

Length  non-dimensionalizing  factor 

16 

GUMULT 

If  1.0  use  U2,  U3  from  continuity  for  right  side 
of  U2",  U3^  equations 

17 

OSG 

Add  stresses  to  RHS  in  pp  equation 

18 

VLDMLT 

Use  laminar  diffusion  in  Uj>  >  U§  equation 

VCMULT 

Add  convection  to  Uj?,  U3  equations 

RARRAY ( 385 ) 

Add  convection  to  UJ  equations 

19 

RHSCAL 

Add  non-symmetric  terms  to  Reynolds  stresses 

U2STRS 

Add  stress  terms  to  right  side  of  pp  equations 

RHOIM 

Adds  wall  damping  to  k,  e 

20 

C4EDSW 

Integration  station  where  ITKE  flag  is  turned  on 

OSUSQ 

TKE  minimum 

21 

SIMPLT 

Station  to  begin  implicit  integration 

TMULT 

Implicit  step  size  =  v  TMULT 

22 

TO 

Initial  integration  station 

TD 

Downstream  distance  to  final  station 

DELP 

Percent  of  TD  for  print 

HSINIT 

Initial  integration  step  size 

23 

CHIEPS 

Convergence  test  on  Newton  Iteration 

HMAX 

Maximum  step  size  allowed 

24 

SEND 

End  of  Namelist  real  data 

50 


III.  DYNAMIC  ARRAY  DIMENSIONING  AND  DEPENDENT  VARIABLE  SPECIFICATION 


25FEDIMN  T  ALLOCATE  STORAGE  FOR  ARRAYS 

24IPINT  -1  T  IiEF'ENBEMT  VARIABLE  STORAGE  ALLOCATION 

27  123547890  0»  000  7*0>  10*11  1  T 


Line  Command 

25  FEDIMN 

26  IPINT 


Dimension  arrays  to  fit  problem  size 

Cards  following  specify  dependent  variable 
and  parameter  arrays 

Variable  1  primary  flow  velocity 

2  Secondary  flow  velocity  (vertical) 

3  Secondary  flow  velocity  (horizontal) 

5  Turbulence  kinetic  energy 

6  Dissipation  function 

7  Perturbation  pressure  (pp) 

8  Complementary  pressure  (pc) 

9  Continuity  equation  potential  (<f>i) 


Dependent  variables  and  parameter  arrays  are 
stored  in  sequential  arrays  each  of  which  is 
NODE  long 


In  NAME01 

NDP  is  the  number  of  arrays  for  space  allocation 

NEQ  is  the  number  of  equations  to  be  integrated 

NEQADD  is  the  number  of  equations  to  begin  inte¬ 
grating  following  initialization 
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IV  FLOWFIELD  GEOMETRY,  Non-Dimensional ization  and  Finite  Element  Matricies 


28LINK2 

14 

T 

GEOMETRY  (DSCRTZ) 

29VX1SCL 

T 

rtSCISSA  COORDINATES 

30  0.0,  4 

.1  .8, 

6 

.6  1.2S,  2  1.  1.  T 

31VX2SCL 

T 

ORDINATE  COORDINATES 

32  0.0,  2 

.4  1.0,  6  . 

9  . 

8,  4  1.  1.25,  4  1.1  .8,  6  1.6  1.25,  2  2.0 

33NDECRD 

T 

DISCRETIZATION  (XI, X 2) 

34  1  25, 

1  13  0 

T 

35ELEM 

0  0  13 

T 

ELEMENT  CONNECTION  TABLE 

3AD0NE 

T 

END  GEOMETRY 

37LINK3 

4 

T 

COMPUTE  NON-DIM.  CONSTANTS.  < DIMEN) 

38LINK1 

3 

T 

FINITE  ELEMENT  MATRICIES  (GEOMFL) 

Line 

Command 

28 

LINK2  14 

Reads  discretization  data  and  generates  triangular 
Finite  Elements 

29,30 

VX1SCL 

Abcissa  Coordinates 

31,32 

VS2SCL 

Ordinate  Coordinates 

33,34 

NDECRD 

Number  of  nodes  in  the  Abcissa  direction 

Number  of  nodes  in  the  Ordinate  direction 

35 

ELEM 

Element  connection  table 

36 

DONE 

End  discretization 

37 

LINK3  4 

Form  non-dimensional izing  parameters 

38 

LINK1  3 

Computes  Finite  Element  matrices 
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V.  TITLES  AND  HEADINGS 


39COMTITLE 

T  READ  TITLE  FOR  COMOC  SYMBOL  PAGE. 

40 DONE 

41DE SCRIPT  204 

T  TITLE  FOR  EACH  OUTPUT-HEADER 

PAGE 

42 

43  CIRCULAR 

JET 

IN  CROSSFLOW. - 3DPNS 

44 

4 5 DONE 

44DESCRIPT  332 

T  PARAMETER  TITLES  FOR  OUTPUT. 

<1  OP  Aft) 

47  REFERENCE 

ENGLISH-FT  ENGLISH-IN 

M~K~S 

C 

48  LENGTH 

...FEET...  ....IN.... 

.  *  # « M . 

.  .  . 

49  VELOCITY 

. .FT/SEC. . 

.  . .M/S. . . . 

•  .  • 

50  DENSITY 

.LBM/FT3. . 

. .KG/M3. . . 

•  .  . 

51  TEMPERATURE 

.RANKINE • • 

.  .KELVIN. . 

52  ENTHALPY 

» BTU/LbM • . 

.  .KJ/KG.  .  . 

53  FROZ.  SPEC. HEAT 

♦ BTU/LBh-R 

.KJ/KG-K. . 

54  VISCOSITY 

.LBN/FT-S. 

.WT-S/ML**  . 

,  .P 

55  LOCAL  PRESSURE 

•  . . PSF  » . • «  ...PSI.... 

. .NT/M2. . . 

.  . . 

54  . . .CPU  TIME. . . 

•MACH.  NO.  ..DPDX1... 

.  .ENERGY. . 

CON 

57  NUGEQM  H'S 

. • , H21 • . . •  ...022.... 

. . . G23* .  .  . 

.  .  • 

58  Xl/LREF 

. DX1/LREF .  .EPSILON.. 

• DX 1H/LREF 

REFL 

REY 

59D0NE 

40MPARA  -1 

T  SCALAR  MULTIPLIERS  FOR  DIMENSIONING  IONUMB 

61 

T  (LOCATIONS  IN  RARRAY) 

62 

5*2’  2*2  142  144  143.  3*2  144 

143,  3*2  170 

174. 

63 

3*2  145  2.  2  -175  3*2.  3*2  174 

2,  3*2  177 

178. 

64 

2  2  149  148  147.  3*2  108 

2*  5*2 p  5*2 

T 

4SI0NUMB  -1 

T  PARAMETERS  TO  PRINT  UNDER  TITLES 

66 

T  (LOCATIONS  IN  RARRAY) 

47 

999.  5*200.  999.  200  4*43.  200 

27  200  2*27 

66 

200  10  200  2*10.  200  58  200  58 

200 1 

69 

200  97  200  97  200.  200  30  200 

30  200 t 

70 

200  38  200  2*38.  999.  39 

4*36 1 

71 

300  154  100  135  122. 

72 

200  184  187  IBS  379,  11  12  14  85  47  T 

73DESCRIP  T  203 

T  TITLES  FOR  PRINTED  ARRAYS  ( IFMTHD > 

74U1  /  UREF 

U2 

/  UREF  U3  /  UREF  PHI 

NU  / 

NUR 

EF  *  RE 

75TKE  /TKEREF 

EPS 

/  EPSREF  PP  /  PSTAG  U1 

PRIME 

PC 

74PPRES 

DPDX 

7  7  DONE 

70IOSAVE  -1  T  ARRAYS  TO  BE  PRINTED.  (CODED  ADDRESSES). 

79  1248  2240  3248  9240  1247 ,  0240  4240  7248  3301  B24B  T 

90 I OMUL T  -1  T  SCALAR  MULTIPLIERS  FOR  EACH  ARRAY  PRINTED. 

81  4*2  21  5*2’  10*1  T 

Line  Command 


39 

COMTITLE 

Problem  identifying  titles 

41 

DESCRIPT 

204 

Std  print  titles 

46 

DESCRIPT 

332 

Std  output  header  labels 

60 

MPARA 

Scalar  parameter  print  multipliers 

65 

IONUMB 

Scalar  parameter  print  locations 

73 

DESCRIPT 

203 

Dependent  variable  print  labels 

78 

IOSAVE 

Dependent  variable  print  locations* 

80 

IOMULT 

Dependent  variable  print  multipliers 

(RARRAY  LOC. 


*Note  in  line  79  that  the  dependent  variables  are  identified  by  decoding  each 
number  specified  (ie.  1248  is  stored  in  the  first  NODE  locations  beginning  at 
the  address  in  IZ ( 248 ) ;  Variable  2,  2248,  is  directly  behind  Variable  1  since 
it  is  stored  in  the  second  NODE  locations  beginning  at  the  address  stored  in 
I Z ( 248 ) ;  etc.). 
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VI.  DEPENDENT  VARIABLE  AND  PARAMETER  BOUNDARY  RELATIONS* 


82KBNO 

1 

83T0P 

-RIGHT 

2 

DONE 

84K8NO 

2 

85T0P 

DONE 

86KBNO 

3 

87LEFT 

TOP 

DONE 

88 KBNO 

5 

89ADD 

DONE 

90 

325*11  1 

T 

91KBN0 

6 

92ADD 

DONE 

93 

325*11  1 

T 

94KBN0 

7 

95T0P 

-RIGHT 

2 

-BOTTOH  2 

DONE 

96KBN0 

9 

97T0P 

-RIGHT 

2 

-BOTTOM  ,  2 

DONE 

Line  Command 


82 

KBNO 

1 

Nodes  where  Ui  is  held  constant  at  initial  values 
(Wall  Nodes) 

84 

KBNO 

2 

Nodes  where  U2  is  fixed  at  initial  values 

86 

KBNO 

3 

Nodes  where  U3  is  fixed  at  initial  values 

88 

KBNO 

5 

Turbulence  kinetic  energy  fixed  at  initial  state 

91 

KBNO 

6 

Dissipation  function  fixed  at  initial  state 

94 

KBNO 

7 

Perturbation  pressure  fixed  at  initial  state 

96 

KBNO 

9 

Continuity  potential  function  fixed  at  initial  state 

Named-side  cards  are  formatted  [4(A8,  18,  4X)] 

Boundary  conditions  are  applied  to  secondary  flow  plane. 


Positive  direction  is  counter¬ 
clockwise  (ie.  [-RIGHT  2]  fixes 
the  right  column  and  node  numbers. 
Node  numbers  are  stored  beginning 
at  the  second  row  from  the  top  and 
proceeding  toward  the  bottom  row). 


X2| 


top 


bottom 


right 

- 

*3 
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VII.  PRIMARY  FLOW  TABLE  DATA 


98VU2P0S  -3 


99 

0. 

0.0166667 

0.166667 

T 

100VU2VAL 

136 

T 

NORMALIZED 

G2.G1 

101 

102VU3P0S 

.5 

-3 

.5  .5 

T 

103 

0. 

0.0166667 

0.166667 

T 

104UU3VAL 

136 

T 

NORMALIZED 

F2,F1 

105 

1. 

1.  1. 

T 

107 

108VPVSX 

109 


43 


PRESSURE 

21*10.05 

PRESSURE 

21*2116.8 


COORDINATES 


COORDINATES 

TABLE  DOWNSTREAM  STATIONS 
0.0  T 

TABLE  CPSF). 

T 


Line  Command 


98 

VU2P0S 

Primary  flow  direction  coordinates  (X^)  where  Xo 
secondary  flow  plane  geometry  transformation  is  defined 

100 

VU2VAL 

Expansion  factors  for  variable  X2  coordinates 

102 

VU3P0S 

Primary  flow  direction  coordinates  (X3)  where  X3 
secondary  flow  plane  geometry  transformation  is  defined 

104 

VU3VAL 

Expansion  factor  for  variable  X3  coordinates 

106 

VX3ST 

Complementary  pressure  data  stations 

108 

VPVSX 

Complementary  pressure  (PSF) 

Geometry  Transformation 

f  Xl 

x2  -  f2l  (Xj) 

n.  =  •  - - — 

£f  2  2  (  X  1  )  ~  fr2l(xx)]/f2 

x  3  ~  f  31  (Xl  ) _ 

[f  32  (Xl  )  -  f3l(Xi)]/f3 

. 


1*2 2  { x i )  and  f  3 2  ( x i )  are  the  span  of  the  secondary 
plane 

f 2 1  (xj  and  f 3 1 ( x i )  are  the  minimum  secondary 
plane  coordinates 

X2  and  x3  are  the  coordinates  of  each  grid  point  in 
respective  directions 

f 2  and  f3  are  the  secondary  plane  normalizing 
lengths 
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VIII.  SOLUTION  INITIATION 


UOLINKCALL 

-1 

T  SUBROUTINES  CALLED  AT  END  OF  GKNUIN 

Ul 

5  12  T  PUNCH  DATA  HEN  LPUNCH  .GT. 

1 12LINK1 

2 

T  GLOBAL  NODAL-ELEMENT  CONNECTIONS 

113RARRAY 

39 

1.0 

114VYY 

115 

78*.  2 

7*. 4  6*. 2  6*. 8  .4  6*. 2 

116 

9  < 

5*1.  .8  .4  6*. 2 

117 

6*.  8 

.4  6*. 2  7*. 4  6*. 2  78*. 2  T 

1 18VYYEND 

1 

119VYY 

120 

325*0.005  T 

121 VYYEND 

5 

122VYY 

- 

123 

325*0.012  T 

124VYYEND 

6 

125R ARRAY 

95 

0.0012 

126R ARRAY 

271 

0.0005 

127READ  5  0 

56  660  T  INITIALIZE  OMEGA  ARRAY 

128 

325*1.0  T 

129LINK1 

4 

T  INITIALIZE  PEP.  VAR. 

130LINK5 

6 

T  SET  NU  TURBULENT  LEVELS.  (DFCFBL) 

1 31IARRAY 

97 

0  T  RESET  ITKE 

1320KNTNT 

T  BEGIN  INTEGRATION. 

133EXIT 

T  TERMINATE  RUN. 

134CASE  END 

T  LAST  DATA  CASE 

Line 

Command 

112 

LINK  1 

2 

Form  global  nodal -element  connection  array 

113 

RARRAY 

39 

Non-dimensional  global  pressure  level 

114 

VYY  1 

Initial  primary  velocity  distribution 

119 

VYY  5 

Initial  turbulence  kinetic  energy  distribution 

122 

VYY  6 

Initial  dissipation  function  distribution 

125 

RARRAY 

95 

Minimum  value  allowed  for  dissipation  function 

126 

RARRAY 

271 

Minimum  value  allowed  for  turbulence  kinetic 
energy 

127 

READ 

56 

Initialize  the  wall  damping  a^ray 

129 

LINK  1 

4 

Initialize  the  dependent  variables 

130 

LINK  5 

6 

Initialize  NU  turbulent  levels 

131 

IARRAY 

97 

Reset  ITKE  to  zero 

132 

QKNINT 

Initiate  implicit  Newton  integration 

133 

EXIT 

Terminate  run  and  return  for  next  data  case 

134 

CASE  END 

Last  data  case 

56 


-C-O-N-O-C 


10IC 
102C 

103  COMMON  /  NPINTQ  /  JJJ(30) 

104C 

105  COMMON  /  VARBLE  /  I ARRAY ( 00500 > f  RARRAY ( 00500 ) 

106  EQUIVALENCE  (  I ARRAY < 00092 > f  IZSIZE  > 

107  COMMON  /  ARRAYS  /  IZ< 130000) 

106  NZSIZE  *  130000 

109C 

HOC  —  IBM  ERROR  HANDLER 

lltc  CALL  ERRSET  <  207f  500*  -1*  If  0*  217  ) 

112  100  CONTINUE 

113  CALL  ZEROTK 

114  CALC  RESET  (  00500 f  IARRAYf  0  ) 

115  CALL  RESET  <  OOSOOf  RARRAY f  0.0  ) 

116  IZSIZE  *  NZSIZE 

117  CALL  RESET  <  IZSIZE*  IZf  0  ) 

113  CALL  DDINPT 

119  GO  TO  100 

120  END 

121  SUBROUTINE  GETPPR 

122  COMMON  /  ARRAYS  /  I Z  < 1 > 

123  DIMENSION  RZ(1>»  L<400) 

124  EQUIVALENCE  <  RZ(l)r  2Z<1)»  L(l)  ) 

125  EQUIVALENCE  < L< 026 > f I INODE > t < L < 109 > f IB1 12  ) r ( L( 110) * IB1 13  ) 

126  EQUIVALENCE  < L ( 077 ) f I AREA  ) t <L<036> r IIELS  ) 

127  EQUIVALENCE  <L  <  089 ) » I XI COR >  #  <L  <  090) * IX2C0R  > 

128  COMMON  /  NPINTQ  /  JUf  JV»  JWf  JHr  JKt  JEf  JP#  JSf  JT#  JAD<21> 

129  COMMON  /  JADRES  /  JUAD,  JVADf  JUADf  JHAD »  JKAD »  JEAD*  JPADf  JSAD 

130  1  *  JEXT ( 22  > 

131  COMMON  /VARBLE/  IARRAY<500)*  RARRAY <500) 

132  EQUIVALENCE  <  I ARRAY < 00014 ) t  NELEH  ) 

133  EQUIVALENCE  <  IARRAY < 00016 > *  NNODE  ) 

134  EQUIVALENCE  <  1  ARRAY ( 00191 > ,  NM  ) 

135  EQUIVALENCE  <  RARRAY < 00027 > f  UINF  ) 

136  EQUIVALENCE  <  RARRAY < 00314 > i  C4  ) 

13?  DATA  NT  /  0  / 

1  38C 

13?  IF  <  NT  .GE.  2  >  RETURN 

140  UNON  *  0.10 

141  UMAX  =  0.4 

142  UMIN  =0.20 

143  YB  *  RZ  < IX2C0R+156  > 

144  ZB  »  RZ<IX2C0R4104> 

145  A  -  1,01  *  <  YB  -  ZB  ) 

146  ASQ  *  A  *  A 

147C 

148C  —  INITIALIZE  DEPENDENT  VARIABLE  ARRAYS 
149C  —  (NON-DIMENSIONAL) 

150C 

151C  JUAD  *  U1 

152C  JVAD  U2 

153C  JUAD  -  U3 

154C  JKAD  «  TKE 

155C  JEAD  ■  DISSIPATION 

156C 

157  DO  200  I  «  I*  NNODE 

158  IM  -  I  -  1 

159  Y  -  RZ< IX2C0R+IH  >  -  YB 

160  Z  -  RZ<  IX1C0R+IM) 

161  rso  *  y*r*z*z 

162  RZ<  JSAD+IM)  ■  1.0 

163  IF  <  RSO  .LT.  ASO  )  GO  TO  200 

164  R  -  SORT  <RSG> 

165  THETA  *  ASIN  < Y/R  > 

166  C0S2TH  -  COS  <  2 • OtTHETA ) 

167  SIN2TH  -  SIN  <  2 . 0* THETA ) 

168  ASR  »  ASO  /  RSO 

169  UVEL  -  ASR  *  ASR 

170  IF  <  UVEL  .LT.  UMAX  )  UVEL  »  UMIN 

171  RZ ( JUAD+IM)  *  UVEL 

172  RZ <  JVAD+ IM )  «  -UNON  *  <  ASR*C0S2TH  +1.0) 

173  RZ ( JUAD+IM )  •  UNON  t  SIN2TH  *  ASR 

174  RZ ( JKAO+ IM )  •  0.1  *  RZ< JKAD+IN) 

175  RZ < JEAD+IM )  ■  0.1  *  RZ(JEAD+IM> 

176  200  CONTINUE 

177  RETURN 

178  END 
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j  79  SUBROUTINE  NOOPPR 

180  COMMON  /  JADRES  /  JUAD*  JVAD*  JWAD*  JHADr  JKAD*  JEAD*  JPAD*  JSAD 

181  l  *  JEXT ( 22 ) 

1 82C 

183  COMNON  /  ARRAYS  /  I Z ( 1 > 

184  DIMENSION  R2(l)f  L<400> 

285  EQUIVALENCE  <  12(1  >*  RZ<1>*  L<1>  ) 

186  EQUIVALENCE  <L< 105) * IPRGRD > 

187C 

188  COMMON  /  VAR6LE  /  I ARRAY ( 00500 > *  RARRAY < 00500 > 

189C 

190  EQUIVALENCE  <  I ARRAY < 000 1 6 ) *  NNODE  > 

191  EQUIVALENCE  (  I ARRAY < 00086 > #  KPNT  ) 

192  EQUIVALENCE  (  I ARRAY ( 00304  >  »  NELPAS  ) 

193  EQUIVALENCE  (  I ARRAY ( 00312 >  »  IUONLY  ) 

19 4  EQUIVALENCE  (  I ARRAY < 00396 > *  JNCPRG  ) 

1  95C 

196  EQUIVALENCE  (  RARRAY ( 00023 ) *  TIME  ) 

197  EQUIVALENCE  <  RARRAY < 00062 ) ►  CONV  ) 

198  EQUIVALENCE  <  RARRAY <001 00 ) f  PPRIME  > 

190  EQUIVALENCE  (  RARRAY < 00207 ) »  CONVRG  ) 

200  EQUIVALENCE  <  RARRAY < 00208 ) »  ANSD2F  ) 

201  EQUIVALENCE  <  RARRAY  <  00347 ) *  TSCALE  ) 

202  EQUIVALENCE  <  RARRAY ( 00385 > *  UCMULT  ) 

203  EQUIVALENCE  <  RARRAY < 00398 ) *  T2PFIX  > 

204C 

205  9600  rORMAT  <  1H  f  6H  PASS  *  14*  1H*  13*  12H  ITERATIONS.  * 

206  1  9H  STATION  *  1PE13.5  ) 

207  c 

..’08  IUONLY  -  3 

209  NVERG  -  CON^G 

210  STAT  «  TIME  *  CONV 

211  NPASP1  *  NELPAS 

212  WRITE  (  6.  9600  )  NPASP1,  NVERG  *  STAT 

213C 

214C  --  GRID  EXPANSION  < CONTRACT  ION > 

215C 

216  CALL  NWGEOM 

217C 

2*18  IF  <  JNCPRG  .EQ.  1  )  CALL  JNCPPR 

219  AMSDIF  -  PPRIME 

220C 

22 1C  —  INCREASE  PRESSURE  EFFECT  TO  100%  OVER  TEN  STEPS 

222C 

223  IF  <  NELPAS  .LT.  9  )  UCMULT  -  1.0E-10 

224  IF  <  NELPAS  .LT.  9  )  GO  TO  100 

225  UCMULT  -  UCMULT  +0.1 

226  IF  (  UCMULT  .GT.  1.0  )  UCMULT  *  1.0 

227  100  CONTINUE 

228  DO  200  I  »  1*  NNODE 

229  IM  »  I  -  1 

230  RZ ( JSA0+ IM )  -  1.0  -  0.5  *  <  RZ < JWAD+IM >**2  +  RZ < JVAD+IH ) **2  > 

231  R2< IPRGRD+IM)  -  UCMULT  *  RZ < IPRGRD+ IM ) 

232  200  CONTINUE 

233  RETURN 

234  END 
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